Draft version February 25, 2013 

Preprint typeset using I^TeJC style emulateapj v. 5/2/11 



THREE-DIMENSIONAL AMR SIMULATIONS OF LONG-DURATION GAMMA-RAY BURST JETS INSIDE 

MASSIVE PROGENITOR STARS 

D. Lopez-Camara^' *, Brian J. Morsony^, Mitchell C. Begelman^'^, Davide Lazzati^ 

Draft version February 25, 2013 

ABSTRACT 

We present the results of special relativistic, adaptive mesh refinement, 3D simulations of gamma- 
ray burst jets expanding inside a realistic stellar progenitor. Our simulations confirm that relativistic 
jets can propagate and break out of the progenitor star while remaining relativistic. This result is 
independent of the resolution, even though the amount of turbulence and variability observed in the 
simulations is greater at higher resolutions. We find that the propagation of the jet head inside the 
progenitor star is slightly faster in 3D simulations compared to 2D ones at the same resolution. This 
behavior seems to be due to the fact that the jet head in 3D simulations can wobble around the jet 
axis, finding the spot of least resistance to proceed. Most of the average jet properties, such as density, 
pressure, and Lorentz factor, are only marginally affected by the dimensionality of the simulations 
and therefore results from 2D simulations can be considered reliable. 
Subject headings: gamma-ray bursts: general — hydrodynamics — supernovae: general 



L INTRODUCTION 

Long-duration gamma-ray bursts (GRBs) are pro- 
duced by collimated relativistic outflows (Sari et al. 1999) 
ejected in the core of massive stars at the end of their 
evolution (Woosley 1993; Hjorth et al. 2003; Stanek et 
al. 2003; Woosley & Bloom 2006). Since their relativis- 
tic outflows have to propagate through their progeni- 
tor star material and exit the star before producing the 
gamma-ray photons, an outstanding issue with this sce- 
nario is to understand the mechanisms that prevent the 
entrainment of baryons in the light, hot jet (MacFadyen 
& Woosley 1999; Aloy et al. 2000). 

On the other hand, even if the jet-star interaction can- 
not slow down the jet, it has a strong impact on its dy- 
namics (Morsony et al. 2007) and can supply enough en- 
ergy to explode the star as a supernova (Khokhlov et al. 
1999; MacFadyen et al. 2001; Wheeler et al. 2002; Maeda 
& Nomoto 2003; Lazzati et al. 2012). In most cases, the 
study of the jet-star interaction has been performed nu- 
merically, with analytic models used only for guidance 
(Aloy et al. 2002; Gomez & Hardee 2004; Morsony et al. 
2007; Matzner 2003; Bromberg et al. 2011). Even so, 
studying the propagation of a relativistic outflow that is 
continuously shocked by a much denser environment is 
not trivial since the length-scale of features in the rela- 
tivistic material is typically ^ R/V and therefore a large 
dynamical range is involved. When possible, adaptive 
mesh refinement (AMR) codes have been adopted (Mor- 
sony et al. 2007, 2010; Lazzati et al. 2009, 2010, 2011b; 
Nagakura et al. 2011), and the simulations have been 
limited to two dimensions (MacFadyen & Woosley 1999; 
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Aloy et al. 2000; MacFadyen et al. 2001; Zhang et al. 
2003; Mizuta et al. 2006; Morsony et al. 2007, 2010; Laz- 
zati et al. 2009, 2010, 2011b; Mizuta & Aloy 2009; Na- 
gakura et al. 2011). These studies have shown that even 
though the jet material is relativistic, the jet-head prop- 
agates sub-relativistically inside the star, thereby allow- 
ing causal contact between the bow shock at the head of 
the jet and the star. The shocked star material there- 
fore drains at the sides of the jet producing a hot cocoon 
(Ramirez-Ruiz et al. 2002; Lazzati & Begelman 2005) 
instead of being entrained in the jet. 

Two dimensional (2D) simulations can provide impor- 
tant answers to the outstanding questions listed above. 
However, they are plagued by artifacts due to the pres- 
ence of a symmetry axis in the center of the jet. First, 
a plug of dense material accumulates in front of the jet 
head, slowing down its propagation and creating plumes 
of hot plasma at wide angles (see Figure 1 in Lazzati et 
al. (2010) for an example). Second, recollimation shocks 
coming from the sides of the jet bounce strongly off the 
jet axis in 2D simulations, while they could dissipate 
more efficiently in a simulation at the natural dimen- 
sionality. Finally, the role of turbulence and instabilities 
cannot be properly explored in 2D simulations. Wang et 
al. (2008) found that in some cases a three dimensional 
(3D) relativistic jet would break apart and not be able 
produce a successful GRB (while in 2D it would produce 
a successful GRB). 

While 3D simulations of GRB jets have been attempted 
in the past (Zhang et al. 2004), they were performed with 
a fixed grid code, casting doubt on their capability to 
resolve the required small scales. A 3D test-case with 
AMR was presented by Wang et al. (2008), but since 
the jet-progenitor evolution varied drastically as a func- 
tion of the numerical resolution (unlike our study), not 
much could be inferred from their study. Thus, in this 
paper we present, for the first time, 3D adaptive mesh 
refinement (AMR) simulations of GRB jets crossing a 
pre-SN progenitor and then flowing through the inter- 
stellar medium. 

This paper is organized as follows. We first describe 
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the physics, initial setup, and the numerical simulations 
in Section 2, followed by our results and discussion in 
Section 3. Conclusions are given in Section 4. 

2. PHYSICS, INITIAL SETUP AND SIMULATIONS 

2.1. Physics and initial setup 

As what now seems to be the generic model used for 
long GRBs (Morsony et al. 2007, 2010; Lazzati et al. 
2009, 2011a,b, 2012; Lopez-Camara et al. 2009, 2010; 
Lindner et al. 2010, 2012; Nagakura et al. 2011, for 
example), we consider the one-dimensional (ID) pre- 
supernova 16TI model from Woosley & Heger (2006) as 
our initial stellar configuration. Initially (in the zero-age 
main sequence) model 16TI is a 16 Wolf-Rayet star 
with 0.01 metallicity, and 3.3 x 10^^ erg s equatorial 
angular momentum. The final outcome of such model is 
a pre-SN progenitor with 13.95 and nearly half the 
size of the sun [Rq = 4.1 x 10^^ cm). Assuming spheri- 
cal symmetry, the ID density and pressure profiles were 
mapped onto a 3D configuration that we assumed to be 
initially without rotation. The internal energy and the 
temperature were calculated assuming a relativistic poly- 
tropic equation of state (7=4/3). The pre-SN progenitor 
was immersed in an interstellar medium (ISM) with con- 
stant density (pism = 10~^^ g cm~^). Even though a 
wind environment would probably be more appropriate, 
we note that within the size of our simulation the dy- 
namical role of the ambient medium is negligible and the 
results are therefore insensitive to the chosen ambient 
medium profile. 

A relativistic jet commencing its flow at the center 
of the pre-SN progenitor was imposed at all times as 
a boundary inflow condition. The jet was launched at 
the center of the star (in fact slightly above it), flowing 
upwards in the polar direction (x=z=0, y=Ri =10^cm). 
The imposed jet had a half-opening angle of 6^=1^^^ a 
constant luminosity of Lq = 5.33x10^^ erg s~-^, an initial 
Lorentz Factor of ro=5, and a ratio of internal over rest- 
mass energy equal to 7^o=80 (Morsony et al. 2007, 2010; 
Lazzati et al. 2009). In order to break the 2D axis sym- 
metry, the jet was slightly asymmetric. For the latter, we 
set the jet with a 1% density and pressure asymmetry on 
either side of a line in the XZ plane 40 degrees from the 
X axis. Diff"erently from Wang et al. (2008) (3D numeri- 
cal study in which a two dimensional symmetrical initial 
setup was assumed) our initial setup resembles that from 
model 3A in Zhang et al. (2004) enhanced with a small 
perturbation in the jet. 

2.2. Numerical simulations 

In order to follow the temporal evolution of our initial 
setup, we solved the 3D gas-dynamic equations using the 
FLASH code (version 2.5) in cartesian coordinates (Fryx- 
ell et al. 2000). The simulation domain covered the top 
half of the pre-SN progenitor star as well as the ISM it is 
immersed in (see for example panel a from Figure 1). The 
boundaries were set at ymin=10^cm, ymax=2.4xl0^^ cm, 

Xmax=-Xmin = 6xl0^^ Cm, and Zmax=-Zmin = 6xl0^° Cm. 

Only the equatorial plane (y=ymin) was set with a reflec- 
tive boundary condition, all the other boundaries were 
set with transmission conditions. We used a 10-level 
binary adaptive grid with square-shaped pixels (Ax = 



A?/ = Az = A). The highest reflnement level (also re- 
ferred to as the flnest resolution level) were accessible 
only at the core of the pre-SN star were the jet is in- 
jected and initially propagates. Moving away from the 
stellar core, the maximum level of refinement was pro- 
gressively decreased. In practice, the base of the jet had 
the finest resolution at all times and the three next finest 
levels followed the jet (and the polar part of the cocoon) 
as it drilled through the progenitor. 

Two set of simulations with a different value of A were 
performed. We will refer to the the hight resolution 
model as "HR", and the low resolution as "LR". The 
HR model had the finest resolution (covering the core of 
the star at all times) equal to A = 3.125 x 10 cm, the 
jet for this case was followed with a resolution of at least 
A = 1.25 X 10^ cm. The LR model had the same setup 
but the value of the the finest resolution was set equal 
to A = 6.25 X 10^ cm, and the jet was followed with a 
resolution of at least A = 2.5 x 10^ cm. The resolution 
with which we follow the jet is comparable to that from 
the 3D collapsar study of Zhang et al. (2004) (where the 
maximum resolution was A ~10^ cm), and to the most 
recent 3D GRB jet study from Wang et al. (2008) (where 
A = 7 X 10^ cm). The resolution with which we resolve 
the core of the star is comparable to that from previ- 
ous 2D GRB jet numerical studies (Zhang et al. 2003; 
Mizuta et al. 2006; Morsony et al. 2007; Nagakura et 
al. 2011). In order to understand the three-dimensional 
effects properly, we also ran an extra two dimensional 
model. 

Differently from the 3D simulation, the 2D run was 
performed in cylindrical coordinates, the polar axis being 
coincident with the jet axis. The 2D model had an initial 
configuration akin to the XY and the ZY planes of the 3D 
model. The 2D model had the same input physics and 
resolution as that of the 3D HR model. A summary with 
the differences between the numerical models is shown in 
Table 1. 



TABLE 1 
Model characteristics 



Model 


A in core 


A in jet 




(xlO'^ cm) 


(xlO^ cm) 


3D LR 


6.250 


2.50 


3D HR 


3.125 


1.25 


2D HR 


3.125 


1.25 



3. RESULTS AND DISCUSSION 

3.1. Global morphology. 

In Figure 1 we show the density stratification maps 
for the XZ, XY, and ZY planes for the 3D LR model. 
Each panel shows a different timeframe: a. t^ = 2.7 s; 
b. th = 4.2 s; c. tc = 5.3 s; d. id = 7.3 s; and e. 
te = 9.3 s. These panels are arranged to illustrate the jet- 
progenitor-ISM temporal evolution (animations of the 
density stratification map, Lorentz factor, radial den- 
sity, radial Lorentz factor, and Schlieren map in the XY 
plane, are linked to the online version of this manuscript). 
The morphology of our system is divided into two main 
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phases: when the jet moving inside the progenitor and 
when the jet has broken out of the star and interacts 
with the interstehar medium. Such temporal evolution 
is consistent with what has already been seen in pre- 
vious numerical studies (Zhang et al. 2003; Morsony et 
al. 2007). Superimposed on the density stratification in 
Figure 1, we show the isocontour levels corresponding to 
10-^, 10-^ 1, 10^, and 10"^ (ah in g cm-^), these isocon- 
tour levels are shown in Figure 2. 

The tbo=4.2 s breakout time is similar to (but some- 
what shorter than) that already seen in previous collap- 
sar studies (Zhang et al. 2003, 2004; Morsony et al. 2007, 
2010). Depending on the progenitor that one chooses, 
and the particular characteristics of the jet, it takes 5 to 
10 s to cross the stellar envelope. Compared to power-law 
stellar models, the models from Woosley & Heger (2006) 
are more compact and dense and it takes less time for the 
jet to cross the realistic progenitors (Mizuta et al. 2006). 
Our tbo is very similar to the breakout time computed 
with the analytical model from Bromberg et al. (2011). 
Still, it must be stated that since the jet in our numer- 
ical simulations is launched at an inner radius which is 
at least 10^ times the gravitational radius (Ri ~ lO^R^ 
for a 1.4M0 black hole), the jet from the simulations is 
somewhat wider than that from the analytical model and 
thus it propagates slower. The tbo for our study implies 
that the average propagation velocity of the jet inside the 
star is ~ 0.32c. The jet, composed of low density mate- 
rial, has its initial opening angle reduced by relativist ic 
hydrodynamic collimation effects. 

Once the jet crosses the stellar envelope and breaks out 
of the surface, the cocoon (which surrounds the jet and 
is present since its formation) expands through the ISM 
(Ramirez-Ruiz et al. 2002; Lazzati & Begelman 2005), 
differently from when the jet is drilling through the pro- 
genitor when the cocoon is bound inside the star and 
close to the jet. When the jet breaks out of progeni- 
tor it becomes uncollimated and the cocoon moves out 
in the polar direction (moving parallel to the jet), also 
expanding sideways on top of the stellar surface. Such 
spreading (see panel c, d, and e from Figures 1-2) was 
predicted by the analytic solution from Bromberg et al. 
(2011). By this time not only does the cocoon present 
zones where variability is clearly present, but also the 
jet presents turbulent-like structures. The variability in 
the cocoon is due to the fact that the jet-cocoon system 
is at least five orders of magnitude denser than the sur- 
rounding ISM. Hence, any instability that forms on the 
cocoon's boundary or that travels upwind from the jet 
into the cocoon is not dissipated. Due to the location 
of the outer boundaries, we are not able to follow the 
jet-cocoon-ISM system entirely after approximately 10 s. 
By this time the cocoon has crossed the outer boundaries 
(the jet crosses the top boundary at approximately 13 s). 
Also, it must be noted that as time passes the inner iso- 
contour {p=10^ cm g~^) disappears. This is due to the 
reverse shock, which is expanding and pushing the dense 
material outwards. Such behavior has already been seen 
in the study of Lazzati et al. (2010). 

3.2. Symmetry loss 

To understand when the cylindrical symmetry is bro- 
ken, in Figure 3 we plot the radial density distribution as 
well as the energy density (U, in erg cm~^) for four dif- 



ferent paths which move along a cone of 2^ half-opening 
angle (with its origin set at x=y=z=0) for model 3D 
LR. One of these paths moves radially (R) along the 
"(-hX,-hZ)" quadrant; another moves in the "(+X,-Z)" 
quadrant; another in the "(-X,-hZ)" quadrant, and fi- 
nally a path which moves along the "(-X,-Z)" quadrant. 

Consistently with previous 2D and 3D collapsar simu- 
lations (Zhang et al. 2003, 2004; Mizuta et al. 2006; Na- 
gakura et al. 2011), we see that as the jet drills through 
the stellar envelope a complex shock system forms, char- 
acterized by a forward and a reverse shock at the head 
of the jet and by and a series of conical recollimation 
shocks. The first recollimation shock, visible almost at 
the base of the jet, seems to be static (at R~2xl0^ cm), 
but this is due to the fact that it moves at relativistic 
speed in the rest- frame. Since the study of the shock 
structure is not our goal, we do not focus on the nature 
of these shocks, nor do we need to know where the con- 
tact discontinuity is set. For the sake of our study all we 
need to be able to discern is the stellar and jet material 
that has and has not been shocked. 

Specifically, the regions which we will be addressing 
to in the rest of the discussion will be the shocked (SJ) 
and unshocked (UJ) parts of the jet. The UJ material 
maintains its initial density profile, while the SJ mate- 
rial breaks the symmetry in the 3D numerical simula- 
tions. The density profile can vary up to two orders of 
magnitude for different locations at the same distance 
from the progenitor center; on the other hand, the UJ 
varies less than an order of magnitude. As the jet crosses 
through the progenitor star its density decreases as a 
function of time (see Figure 4). Before the jet breaks 
out from the stellar surface, the density profile inside 
the progenitor follows a quasi- const ant profile which for 
tbo is ~10~^ g cm~^. Then, when the jet breaks out 
of the stellar surface, it recovers a decaying radial den- 
sity profile that for 9.3 s reaches density values as low as 
10"^ g cm"^. 

3.3. Lorentz factor evolution 

In Figure 5 and Figure 6 we show the temporal evo- 
lution for the Lorentz factor (with the velocity field also 
present); and the radial Lorentz factor profile along the 
2^ radial path for model 3D LR. Before the breakout 
time only a relativistic jet (with F ^10) is present (see 
panel with t=4.2 s from Figure 5). In Figure 6 we see 
that the SJ material for t<tbo (blue, red and green lines) 
reaches values close to F =15; and the UJ Lorentz factor 
remains practically the same as the initial Lorentz fac- 
tor (ro=5). This behavior is consistent with what has 
already been seen in previous GRB jet numerical sim- 
ulations where the initial Lorentz factor, prior to tbo, 
reaches values close to 10 (Zhang et al. 2004; Mizuta et 
al. 2006). Once the jet breaks out of the stellar surface, 
the jet is accelerated. The high internal energy is able 
to accelerate material with Lorentz factors values of or- 
der F ^100 in some zones. If accelerated with no energy 
dissipation, the jet's maximum Lorentz factor would be 
400 (Foo=ror?o) (Morsony et al. 2007). In the panel with 
t=5.3 s from Figure 5 (brown line in Figure 6) we show 
how the jet's forward shock and the recently formed co- 
coon produce a "mushroom-like" high-F structure. At 
later times (t>tbo) (orange, cyan and black lines in Fig- 
ure 6) the mushroom-like structure grows bigger and its 



4 



Lopez-Camara et al. 




Density 


(g cm"3) 


1x105- 


■ 


3x102- 




1x100- 




3x10-3- 




1x10-5- 





c. 5.3 s 



Fig. 1. — Density stratification maps (g cm ^) for different timeframes (a. 2.7 s; b. 4.2 s; c. 5.3 s; d. 7.3 s; e. 9.3 s) for model 3D LR. 
The isocontour levels correspond to: 10^ g cm""^; 10^ g cm~^; 1 g cm""^; 10~^ g cm""^; 10~^ g cm~^. In order to better visualize the 
internal structure in the pre-SN, the minimum value in all the density stratification plots was set to 10~^g cm"*^. A movie of this figure is 
available in http:/ /www4.ncsu.edu/~dlopez/Simulations_(published).html 



Lorentz factor F increases significantly. 

To see the high Lorentz factor material in the jet, in 
Figure 7 we plot the F isocontours for t=9.3 s. By this 
time the mushroom structure is evident and also certain 
regions in the polar axis reach Lorentz factor values as 
high as F = 50 (pink isocontours). Unfortunately our nu- 
merical domain does not permit us to follow such high-F 
regions and they escape the top boundary after approx- 
imately 10 seconds. This is once more congruent with 
the results from previous studies where after the tbo the 
cocoon reaches values as high as F ~15, and the jet val- 
ues of order F ~100 (Mizuta et al. 2006). It must also be 
noted that when the jet breaks out of the stellar surface, 
a low- speed wind forms. This wind expands isotropically 
from the point in the stellar surface where the jet drilled 
through, and moves at an average speed vastly inferior 
(v<0.01 c) to that of the jet. 

3.4. Resolution effects 



In order to be able to evolve the initial setup up to 
integration times of order ^10 s and to resolve the jet- 
progenitor with a suitably fine grid an AMR mesh was 
used. In Figure 8 we show the fraction of the volume that 
the three finest resolution levels occupied as a function of 
time. The finest grid level, the one with which the base 
of the jet was resolved (A ~10^ cm, red line in Figure 8) 
occupied less than 10~^ of the entire volume. Meanwhile, 
the two next finest levels, which followed the propagation 
of the jet through the progenitor (A ~10^ cm, blue and 
green lines in Figure 8), occupied less than 10~^ and 
10~^ of the volume. Needless to say, if we had used a 
fixed mesh with a comparable resolution to that with 
which the base of the jet was resolved, it would have 
required ~10^ more computational power (compared the 
computational power used in our simulations), and thus 
the benefit from using an AMR scheme. 

To verify that the evolution of the jet from our results 
is not dependent on the numerical resolution, we ran a 
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Fig. 2. — Density stratification map for different isocontour levels (red=10^ g cm ^; yellow=10^ g cm ^; green=l g cm ^; 
cyan^lO"-^ g cm~^; blue=10~^ g cm~^) for model 3D LR. The timeframes are the same as those indicated in Figure 1. A movie of 
this figure is available in http://www4.ncsu.edu/~dlopez/Simulations_(published).html 



new model with the same setup and physics but with a 
maximum resolution two times finer than for the low res- 
olution (see Section 2.2 for details). In Figure 9 we show 
the density profiles for the 3D HR and 3D LR case. In 
each case we show the two main phases already discussed 
in Section 3.1: prior to the breakout phase; the break- 
out; and the post breakout phase. We must note that the 
selected timeframes for each of the phases were chosen 
arbitrarily so that the LR and HR cases resemble each 
other, and hence their basic morphology characteristics 
can be compared. The main result is that the basic mor- 
phology characteristics from each phase (see discussion 
in Sections 3.1—3.3), are well reproduced independently 
of the numerical resolution. Unlike the results from the 
numerical 3D jet GRB study from Wang et al. (2008) 
where high resolution models gave qualitatively different 
jet dynamics, we obtain consistent jet behavior indepen- 
dently of the resolution. 

Among the differences associated with the resolution, 
are a higher level of turbulence and a slower advance of 



the jet head in the HR model. The latter's jet moves 
~20% slower than the LR case, hence the breakout time 
for the HR case is tbo=5.1 s. The jets velocity resolution 
difference is due to the fact that the HR case has a wider 
jet (~5% wider than the LR case). Since we are power- 
ing both jets equally, the narrow-LR jet will move faster. 
The turbulence resolution difference is due to the fact 
that the LR simulation has higher diffusion, and thus 
suppresses the small scale instabilities which are present 
in the HR model. The higher amount of turbulence in 
the HR model also slows it down (compared to the LR 
model), this due to the fact that a larger fraction of the 
energy is converted into turbulence. Hence, reducing the 
HR jet's kinetic energy and ram-pressure. We must note 
that these two resolution effects are consistent with what 
has already been seen in previous jet-collapsar simula- 
tions, for example Morsony et al. (2007). Even though 
the latter study is a two-dimensional one, it also presents 
more vortices in the HR than in the LR case. 
To illustrate how in the HR model there is more vari- 
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Fig. 3. — Radial profiles for different 2^ paths and times for 
model 3D LR. Panel a through d show different radial density 
profiles (g cm~^), panel e shows different energy density pro- 
files (erg cm~^). Each of the radial paths commences in the 
origin and runs through different quadrants: (+X,+Z) quadrant 
(red line); (+X,-Z) quadrant (green line); (-X,+Z) quadrant (blue 
line); (-X,-Z) quadrant (black line). Each panel corresponds to 
a different timeframe: a. 2.7 s, b. 3.5 s, c. 3.9 s, d. and e. 
4.2 s. The forward and reverse shock Mach number for each time- 
frame are also indicated. A movie of this figure is available in 
http: / / www4.ncsu.edu / ~dlopez / Simulations_(published) .html 



ability than in the LR one, in Figure 10 we plot the 
radial density profiles for both resolution models. The 
upper panel of Figure 10 is the radial density profile for 
a path within the jet (specifically a 2^ path), while the 
lower panel is the radial density profile in the edge of the 
jet (10^ path). Inside the jet, there is no major difference 
between the two resolution models. On the other hand, 
at the edge of the jet the numerical resolution clearly af- 
fects the density profile. Here the HR presents numerous 
depressions in the density radial profile while the LR case 
has a radial profile which follows a smoother distribution 
with less depressions and variability. 

We also analyzed the effects of numerical resolution 
on the jet's Lorentz factor distribution. Before the jet 
breaks out, apart from the velocity with which the jet 
evolves, there is no clear difference between the low and 
high resolution models. But at t>tbo there are morpho- 
logical changes due to the resolution. Figure 11 shows 
the Lorentz factor map for the 3D LR and the 3D HR 
models just after breaking out of the stellar surface. The 
Lorentz factors are comparable but, as expected, the HR 
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Fig. 4. Time evolution (0 s; 5.3 s; 7.3 s; 9.3 s) for the 
2° radial density profile (g cm~^) from the (+X,+Z) quadrant 
(black line in Figure 3). The stellar surface is indicated by 
the cyan dashed line. A movie of this figure is available in 
http: / /www4.ncsu.edu/~dlopez/Simulations_(published).html 



Lorentz 
factor 




6-6 
x(10io cm) 

Fig. 5. — Lorentz factor stratification maps and velocity 
field for different timeframes (4.2 s; 5.3 s; 7.3 s; 9.3 s) 
for model 3D LR. The brown dashed line indicates the 
stellar surface. A movie of this figure is available in 
http: / /www4.ncsu.edu/~dlopez/Simulations_(published).html 



case has more turbulent-like structures. For both resolu- 
tions the low-density material (situated along the polar 
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Fig. 6. — Same as Figure 4 but for the Lorentz 
factor. A movie of this figure is available in 

http: / / www4.ncsu.edu / ~dlopez / Simulations_(published) .html 
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Fig. 8. — Temporal evolution of the fraction of the volume that 
the two finest level occupied. The red line corresponds to the finest 
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Fig. 7. — Lorentz factor isocontour map at t=9.3 s for 
model 3D LR. The isocontour levels correspond to: 2, 5, 
10, 20 and 50. A movie of this figure is available in 
http: / / www4.ncsu.edu / ~dlopez / Simulations_(published) .html 



axis) has higher Lorentz factors compared to the material 
in the edge of the jet. The main difference, in quahtative 
terms, is that the high-F regions from the LR model are 



split into smaller regions with lower F values in the HR 
model. 

Another interesting aspect of our results is that the 
jet Lorentz factor morphology resembles the precessing 
jet case from Zhang et al. (2004). In our simulations 
though, the zig-zagging pattern in the Lorentz factor 
seems to be due to the ability of the high-F material 
to wobble around the star's rotation axis and propagate 
through paths of least resistance (see below for more de- 
tails). The observation of this effect is likely facilitated 
by the fact that we have a larger dynamic range inside 
the star {Sr = Ro/Rmin = 41) compared to the one from 
Zhang et al. (2004) {Sr = 8.8). Thus, there is enough 
range inside the star for 3D instabilities to develop. In 
addition, the individual grid pixels in the Zhang et al. 
(2004) simulations was not square. Rectangular pixels 
generate more diffusion in the longest of the pixels di- 
rection, and results are therefore not as robust as those 
from square pixel grid simulations (where the diffusion is 
the same for all three directions). 

Finally, we must remark that we do not claim to reach 
convergency. If we take the number of grid cells across 
the jet diameter as an estimate of the Reynolds number 
(Re) of the simulation, we see that at the present time 
we can only reach Re~200. Such Reynolds number is ap- 
proximately two orders of magnitude below the required 
Reynolds number in which the jet behavior is indepen- 
dent of the resolution (Birch & Eggers 1972). Unfortu- 
nately, numerical simulations as those presented in this 
study, with resolutions of at least two orders of magni- 
tude finer are not feasible (due to technical difficulties) 
to date. 

3.5. 2D vs 3D simulations 

Finally, we checked how the evolution of the jet 
through the stellar envelope varies in two and three 
dimensional simulations. For this we ran an ex- 
tra 2D numerical model with the same resolution 
(A=1.25xl0^ cm), and the same parameter values (lu- 
minosity, Ri, 6>o, Fo, and r]o, see section 2.2 for more 
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Fig. 9. — Density stratification maps (g cm~^) for the 3D LR model (upper panels) and the 3D HR model (lower panels). For both 
models, we show representative timeframes from each of the two main phases (see text for discussion). The isocontour levels are the same 
as the ones indicated in Figure 1. 
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Fig. 10. — Radial density profile (g cm~^) for the 3D LR model 
(red line) and the 3D HR model (black line). For both resolutions 
the 2^ radial paths from the (+X,+Z) are shown in the upper panel, 
meanwhile the lO'^ paths are shown in the lower panel. 



details) as the 3D HR case. Apart from the dimension 
difference from the 3D models, we assumed polar axis 
symmetry in the 2D simulation; thus in reality we only 
simulated half of the x-axis domain (i.e Xmin=0). 

The 2D simulation was carried out in cylindrical coor- 
dinates, with the polar axis coincident with the jet axis. 
In Figure 12 we show the density stratification maps for 
the 2D model at the two phases (as well as for the break- 
out time): a. t<tbo; b. t^tbo, and c. t>tbo)- The time- 
frames for each of the 2D phases shown in Figure 12 were 
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Fig. 11. — Lorentz factor stratification maps for the 3D LR (left 
panel), and the 3D HR model (right panel). For both resolution 
models, we show a representative timeframe for when the jet has 
already broken out of the star and is moving across the ISM. 



chosen so that they resemble the correspondent time- 
frames from the 3D HR case. The basic morphology 
in the 2D case resembles that from the 3D model. In 
both cases we see a collimated jet that manages to drill 
through the stellar envelope. Apart from the polar axis 
symmetry imposed in the 2D model, there are many sub- 
tle differences between the 2D and 3D results: 
1. The jet moves slower in the 2D model than in the 



3D GRB simulations 
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3D one (congruent with Zhang et al. (2004)). Thus, the 
2D breakout time is larger (t^^=7 s) than for its corre- 
spondent 3D HR model. The reason for the slower jet's 
motion in the 2D model is the imposed symmetry. Not 
only is the jet axis- symmetric, but also the stellar ma- 
terial in front of the jet has to remain symmetric at all 
times. The SJ material can only escape from the jets plug 
sideways, so a lot of energy ends up going into accelerat- 
ing this stellar material. In 3D models, instead, the jet 
can deflect slightly and go around the plug (rather than 
continuing to accelerate it, see below). Finally, the 2D 
jet model has a SJ material mildly broader (?^20%) than 
its respective from the 3D model. 

2. Even though the 2D model has the same resolution, 
the 2D jet presents less turbulent-like morphology than 
what is present in the 3D HR case. Also, the cocoon is 
less turbulent and broader in the 2D scenario, as is the 
case in the numerical study of Zhang et al. (2004). 

3. Two low-density plumes are present in the 2D simula- 
tions (see right panel from Figure 12). It must be noted 
that due to the imposed axis-symmetry, the plumes ac- 
tually correspond to a low-density torus around the jet 
head (if it were a three dimensional domain and not two 
dimensional). Such low-density torus is not present in 
any of the 3D simulations. This is somewhat similar to 
the findings from Zhang et al. (2004) where the head of 
the jet is noticeably different depending on the simula- 
tions dimensionality (either 2D or 3D). 

4. As shown in the upper panel of Figure 14, where we 
plot the density profile along the polar axis for the 2D 
and 3D models, the 2D density radial profile is less tur- 
bulent and the depressions are more profound than those 
in the 3D density radial profile. Also, the SJ material is 
mode dense by nearly two orders of magnitude in the 2D 
scenario close to the jet head. 

In Figure 13 we show the Lorentz factor structure for 
the 2D case (once the jet has just broken out of the stellar 
surface). Even though the mushroom F structure forms, 
the 2D Lorentz factor morphology is noticeably differ- 
ent to that from the 3D HR case. The 2D F structure 
presents much less variability. 

The 2D low-density regions have high F values of up 
to 15-20, values which are in agreement with those ob- 
tained by other 2D GRB jet studies (Zhang et al. 2004; 
Nagakura et al. 2011). The 2D model also has a SJ which 
is broader than that from the 3D case. As was the case 
for the density map in the 2D model, the head front of 
the cocoon has less turbulent-like F structures. In order 
to clarify this point, we show the radial Lorentz factor 
profile (along the polar axis) in Figure 14 (lower panel). 
Not only is a smoother radial profile present in the 2D, 
but also the high-density low-F) relationship is present. 
The SJ from the 2D case is approximately two orders of 
magnitude more dense than from the 3D model, but has 
a rather smaller F value (of at most 5). 

In order to analyze how much the jet changes direction 
as it drills through the progenitor star (and later through 
the ISM once the jet has broken out of the star) in a three 
dimensional domain, we plot in Figure 15 the energy den- 
sity {U) map in the XZ plane. The XZ planes shown for 
each timeframe correspond to the position where the U 
centroid of the forward shock front was located at (see 
caption of Figure 15 for more details). Panels a through 



d show the U map for when the jet is drilling through the 
stellar progenitor. In these it is noticeable how the cen- 
troid of the forward shock (CFS) does not have a gaussian 
like profile (it may have turbulent like behavior or even 
multiple spikes) and how the CFS wobbles around the 
polar axis finding the spot of least resistance to proceed. 
For example, notice how just before the jet breaks out 
of the star (panel d), the CFS is located far from the 
polar axis (x=-0.3xl0^^ cm , z=-0.6xl0^^ cm). Panel e 
shows how once the jet has broken out of the star and 
the cocoon has expanded thoroughly around the progen- 
itor star, its correspondent CFS also expands and also 
remains far from the polar axis. 

To further understand the defiection of the jet inside 
the pre-SN progenitor, we show the temporal evolution 
of the angle between the CFS and the polar axis (^, black 
line in Figure 16). For a jet that is well aligned with the 
polar axis then the CFS displacement angle would yield 
^ = 0, clearly in Figure 16 this is not the case and the 
jet wobbles inside the star (with oscillating between 
0.1^ and 2^). Hence, the jet moves faster in 3D than 
in 2D because it is able to wobble and move along the 
path with least resistance (apart from having a narrower 
jet-cocoon). Note that is always within the relativistic 
collimation angle {I/O, red line in Figure 16), thus the 
relativistic jet is causally connected at all times. 

3.6. Limitations and comparison to other work 

As with all numerical work, the choices made in carry- 
ing out the simulations refiect intentions and biases, and 
the current investigation lacks in several aspects. For 
example, similar to Zhang et al. (2004), Morsony et al. 
(2007, 2010), and Lazzati et al. (2009) we assumed that 
the star was static at all times which is clearly not the 
real case as the pre-SN for long GRBs have very high 
angular momentum values (J>10"'^^ cm^ s) (Woosley & 
Heger 2006). We justify this by pointing out that the 
dynamical timescale of the pre-SN is of order close to 
hours. Then, since the integration time in our numerical 
simulations was of order 10 s, we were safe to assume that 
the pre-SN progenitor remained practically static at all 
times. In a previous study with a similar setup (Lazzati 
et al. 2011b) found that after 10^ s the pre-SN stehar 
envelope had only expanded 2% of its original size. 

Another issue which can be improved is the ISM distri- 
bution. The pre-SN progenitor that we use as the initial 
setup has no hydrogen shell since during its stellar evo- 
lution it was lost by the presence of a stellar wind (which 
will also affect the ISM surrounding the pre-SN star). So, 
in order to have full consistency the ISM should have a 
density profile which was affected by the pre-SN wind, i.e. 
a profile that follows a ocR"^ distribution (Zhang et al. 
2003, 2004; Cannizzo et al. 2004; Nagakura et al. 2011). 
But since the jet-cocoon system is an ultra-relativistic 
fiow, the density profile of the ISM will barely affect the 
jet once this has just broken out of the stellar surface 
(Morsony et al. 2007). In fact, the GRB-jet needs to 
reach ^lO^^cm for the ISM's profile to play a key role 
in the fiow (Blandford & McKee 1976; De Colle et al. 
2012). Thus, we were secure to assume that the ISM 
density was constant. 

We use an adiabatic (7=4/3) as our equation EOS 
(Zhang et al. 2003, 2004; Mizuta et al. 2006; Morsony 
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Fig. 12. — Density stratification maps (g cm ^) for the 2D model. For this we show representative timeframes from each of the two main 
phases (see text for discussion). The isocontour levels are the same as the ones indicated in Figure 1. 
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Fig. 13. — Lorentz factor stratification map for the 2D model at 
t=7.7 s. 



et al. 2007, 2010; Lazzati et al. 2009; Nagakura et al. 
2011). We do not take into account the neutrino pres- 
sure, nor do we take into account the gravitational ef- 
fects from the central compact object. Even though it 
has been shown that close to the pre-SN's progenitor 
nucleus the neutrinos play an important role (Lopez- 
Camara et al. 2009), since the inner boundary was set 
so far away, Ri ^lO^cm, equivalent to approximately 10^ 
gravitational radii, from the region where neutrinos dom- 
inate (and where the compacts object relativistic effects 
must be taken into consideration), the neutrino and rel- 
ativistic effects were safely ignored. 



2 4 
R (10^0 cm) 

Fig. 14. — Radial density profile (g cm""^) (upper panel), and 
the radial Lorentz factor (bottom panel) for both the 2D model 
(red line) and the 3D model (black line) for the timeframe when 
the jet has broken out of the star. For both models the path is a 
polar axis (0^) radial paths from the (+X,+Z) quadrant. 



Since the follow up of newly formed elements was not 
the aim of this study and that the calculation of such 
new elements would have not permitted us to study the 
flow both at an adequate resolutions and for the long 
times desired, nuclear burning was not included. Also, 
even though magnetic fields will affect the emissivity of 
the jet (Mizuta et al. 2006), and could even give ori- 
gin to variability in the light curve (Balbus & Hawley 
1998), they were disregarded due to the technical diffi- 
culties when following a magnetized relativistic flow with 
an adaptive mesh code. 

4. CONCLUSIONS 

We present, for the first time, 3D AMR simulations 
of ORB jets expanding inside a realistic pre-SN progen- 
itor and then flowing through the interstellar medium. 
Our numerical simulations, confirm that relativistic jets 
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Fig. 15. — Energy density (erg cm^) XZ stratification maps for 
the centroid of the head front of the jet-cocoon structure for dif- 
ferent times. Panels a through e correspond to the 3D HR model, 
panel f correspond to the 3D LR model. Panel a. t=ls, b. t=2s, 
c. t=4s, d. t=5.1s, e. t=5.7s, f. t=10s; have the XY plane located 
at Y/(10io cm)=1.3, 3.3, 18.9, 39.9, 57.1 (respectively). In each 
panel the maximum and minimum values of the forward shock's 
energy density centroid (in erg cm^) is indicated. Notice how some 
panels have different scales. 




Fig. 16. — Temporal evolution of the CFS (black line), and the 
relativistic collimation angle (red line). 



can propagate and break out of the progenitor star while 
remaining relativistic. 
The morphology is divided into two main phases: 

1. Pre-t^o- During this phase the jet head moves at 
mildly relativistic velocities (~c/2) inside the progeni- 
tor's stellar envelope. 

2. Post-t^o. Once the jet breaks out of the surface, it 
accelerates and reaches Lorentz factors of order F ~ 50. 

The initial progenitor density profile is reshaped by the 
forward and reverse shocks. The material between the 
forward and reverse shocks break the two-dimensional 
symmetry in the numerical simulations. 

We obtain similar behavior independently of the nu- 
merical resolution. The resolution does not affect in great 
detail the flow and the morphology in each phase is well 
reproduced. Still, the amount of turbulence and variabil- 
ity observed in the simulations is higher for higher resolu- 
tions. Also, for finer numerical resolutions the jet moves 
slower; and regions with high Lorentz factors break up 
into smaller regions with lower F values. 

The propagation of the jet head inside the progeni- 
tor star is slightly faster in 3D simulations compared to 
2D ones at the same resolution. This behavior is due 
to the fact that the jet in 3D simulations is narrower 
and can wobble around the jet axis finding the spot of 
least resistance to proceed. Most of the jet properties, 
such as density, pressure, and Lorentz factor, are only 
marginally affected by the dimensionality of the simula- 
tions and therefore results from 2D simulations can be 
considered reliable. If, instead, more detailed proper- 
ties such as variability are to be investigated, simulations 
carried out in the proper dimensionality (i.e. 3D) are re- 
quired. 
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